Wavevector resonance in a nonlinear multi-wavespeed chaotic billiard 
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Nonlinear coupling between eigenmodes of a system leads to spectral energy redistribution. For 
multi-wavespeed chaotic billiards the average coupling strength can exhibit sharp discontinuities 
as a function of frequency related to wavevector coincidences between constituent waves of dif- 
ferent wavespeeds. The phenomenon is investigated numerically for an ensemble of 2D square 
two-wavespeed billiards with rough boundaries and quadratic nonlinearity representative of elasto- 
dynamic waves. Results of direct numerical simulations are compared with theoretical predictions. 



Recent years have been marked by increased attention 
to multiple scattering and propagation of classical waves 
in disordered nonlinear media. Main research has been 
in optics; see 0] and references therein. Characteristics 
of a wave-field pertinent to continuous- wave problems, 
such as angular correlations and coherent backscatter- 
ing, have dominated [2| . Influence of nonlinearity on the 
stability of a speckle pattern has been considered also Q; 
nonlinear phenomena in transient fields, as, for example, 
frequency shifting of a lasing-mode in a random laser Q, 
have been studied. 

An anisotropic or multi-wavespeed nature of the 
medium supporting propagation of the classical waves al- 
lows a broader ground for interplay between nonlinearity 
and multiple scattering. Statistical effects of nonlinear 
wave-field behavior have received little attention so far. 
A class of systems in which these effects are expected, is 
given by elastodynamic waves in solids, with anisotropy 
or multiple speeds of propagation as a rule rather than 
exception. Ballistic billiards in the form of elastic solids 
with a ray-chaotic shape are conveniently realizable and 
allow readyaccess to the time domain and observation of 
transients Q . These billiards are representative of classi- 
cal wave-bearing systems, and are suited for the study of 
statistical nonlinear effects that stem from the multiple- 
wavespeed character of the wave-field. The purpose of 
this study is to provide evidence of one of such effects, 
namely wavevector resonance in a chaotic nonlinear bil- 
liard, with an elastodynamic billiard chosen as an under- 
lying physical model. 

Dynamics of an elastodynamic billiard after excitation 
by a transient source can be reduced in the absence of 
energy dissipation to a set of nonlinearly coupled oscilla- 
tors with natural frequencies ujj,, amplitudes of vibration 
dk, and coupling matrices N 

d k + uldk + Nkimdidm + N k i mn did rn d n = 0, (1) 

where nonlinearity up to cubic terms has been accounted 
for. Associated with each oscillator (mode) is its lin- 
ear energy E k = [d\ + bj\d\^j /2. Although not being 

true energy in the presence of nonlinearity, E k elucidates 
trends of modal energy redistribution when nonlinearity 
is weak. 

For observation times much smaller than the inverse 



modal spacing, t -C D, individual modes are not resolved, 
and statistical description is in order. One of the energy 
quantities conveniently accessible for experimental mea- 
surement under these conditions is the average spectral 
density, E (w k ,t) — D (cu k ) (E k (t)), a constant in the ab- 
sence of dissipation or nonlinearity. Under the influence 
of nonlinearity it evolves in time. The evolution leads to 
energy deposit into frequencies that contained no energy 
initially, thereby allowing detection of the nonlinearity. 
In case of a weak nonlinearity the leading contribution 
to the deposit is given by convolution of energy densities 
at two frequencies that have a given (target) frequency 
w as a combination, i.e. their sum or difference [f|: 

r+ao 

E{u,t) = D{uj) I duj' VN(w,w', \u±uj'\) 
Jo ± 

xE{\w±u'\)E{w')/(oj±jfu' 2 . (2) 

Redistribution of the energy involving triads of frequen- 
cies is characteristic of a dominant quadratic nonlinear- 
ity - contribution of the cubic terms of eq. (QJ vanishes 
on the average, being of the next order of smallness, - 
- and is in agreement with behavior expected from ele- 
mentary theory of nonlinear oscillations (cf., for example, 
0|) • It corresponds to internal (frequency) resonance be- 
tween the modes of the billiard, which effectively takes 
place when individual modes are not resolved, and hence 
a combination frequency of two source modes is indis- 
tinguishable from the frequency of another mode of the 
system. The resonance manifests itself in early-time lin- 
ear energy growth that can be linked to secular terms 
arising in Regular Perturbation Theory Q. 

Redistribution of the spectral energy in the billiard oc- 
curs due to nonlinear coupling of the modes. The average 
coupling strength is provided by the function, 

N(u> k ,u>i,u} m ) = (7r/2) (N klm N k i m ) 

= (tt/2) N al3l N vm K av (uj k ) Kp^ (uji) K in (u m ) . (3) 

Greek indices denote spatial degrees of freedom includ- 
ing both spatial position x and Cartesian indices, and 
repeated indices imply summation. The function is sym- 
metrical with respect to its arguments, rendering inter- 
action strength of any triad of frequencies independent 
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of energy transfer direction between them. Nevertheless, 
in nondispersive billiards with the inverse modal spac- 
ing given by Weyl series, D ~ VdUJ^ 1 + . . . jg], where 
Vd stands for d-dimensional volume of the billiard, the 
overall power input exhibits a trend of the energy to be 
transfered up the frequency spectrum. 

Specifics of the physics of a particular system re- 
side in statistics of its modes Uk, and operator TV that 
acts on them, Nkim = N a ^Uk (a) ui u m (7)- Un- 
der the assumption that the modes are mutually uncor- 
rected Gaussian random vectors (3, UJJ|, the statistics 
are fully provided by the pairwise spatial modal correla- 
tor K a f3 (ujk) = (tife (a) %(/?)). In the short-wavelength 
limit, when boundary effects are unimportant, the corre- 
lator can be approximated by means of the Green's func- 
tion in boundless medium, with long-range correlation 
being limited to billiard diameter L. Refined approxima- 
tions of the correlator that include finite size andgeom- 
etry effects can be constructed upon need, see lllll and 
references therein. 

Under a Random Wave Model (RWM) a mode in- 
side a chaotic billiard can be viewed as superposition 
of constituent plane waves coming from all directions 
and having random amplitudes and phases jjj. In a 
multi-wavespeed billiard these waves can have different 
speeds of propagation and different amplitudes, requir- 
ing a modification to RWM [12J. Each pair of the con- 
stituent waves from the source modes interacts nonlin- 
early producing a constituent wave of the target mode. 
Two conservation laws must be obeyed in the process. 
The first requires frequency of the resultant wave to be 
equal either to the sum or difference of the source ones. 
A parallel can be drawn to energy conservation applied 
to particle scattering problems. The law is automati- 
cally satisfied by the internal-resonance structure of eq. 
0. The second law requires wavevectors of the con- 
stituent waves to sum, thus standing for conservation of 
wave pseudo-momentum. In nondispersive systems in- 
teraction of constituent waves of the same wavespeed is 
always allowed by the above conservation laws, when the 
waves are collinear. This mechanism provides a non-zero 
background coupling strength for any frequency combi- 
nation of source modes. Interaction of constituent waves 
of different wavespeeds, however, is only possible if the 
frequencies (and polarizations) of the source waves are 
in special relation to each other If the frequencies 

are such, this additional interaction channel is open, and 
nonlinear coupling strength is expected to be greater. In- 
tegral contribution of all constituent- wave interactions is 
provided by eq. QJ, statistics of the waves being in- 
corporated by means of spatial correlator K. Precise 
conditions under which the maximum of the coupling 
strength is achieved depend on the specifics of the non- 
linearity; for typical elastic solids they were seen to im- 
ply interaction of the constituent plane waves in a nearly 
collinear fashion By analogy to internal (frequency) 
resonance a sharp increase in the coupling strength due 
to nonlinear interaction of constituent waves with differ- 



ent wavespeeds is termed wavevector resonance here. 

To verify existence of the wavevector resonance a series 
of direct numerical simulations (DNS) was performed. 
An ensemble of discrete 2D square billiards with aver- 
age boundary roughness of one twenty-fourth of the bil- 
liard size was taken; cf. 0|- Evolution of the wave- 
field inside the billiards was governed by first-order finite- 
difference version of the following two-wavespeed elasto- 
dynamic equation: 

Ui = CfUtjj + (cf - Cf ) Uj ti j + S$ijklmn (uk,lU m ,n) j , (4) 

where cj and ct are longitudinal and transverse 
wavespeeds, respectively. The nonlinear coupling term 
was chosen to be quadratic in derivatives of the displace- 
ments u, a form representative of the physical nonlin- 
earity in isotropic elastic solids described by the five- 
constant theory ^^|. The derivatives were coupled by 
the elementary isotropic tensor, 

^ijklmn — 

g(fiik5jmSln + Sik5jnSl m + 5u5j m 5kn + 5u5j n 5km 
+5i m 5jk5ln + SimSjlSkn + 5i n Sjk5i m + S in Sjl5km)- 

Strength of the nonlinearity was controlled by param- 
eter e, which was kept small on the order 10~ 2 , while 
u = 0(1). Dirichlet boundary conditions were imposed 
on the billiard boundaries. Initial displacement field was 
arranged so that only modes supporting two narrow- 
band Gaussian peaks of the spectral energy centered at 
given frequencies u>i and u>2 were excited. The width of 
the peaks was taken wide enough to encompass tens of 
modes. Evolution of the initial field up to one tenth of 
the Heisenberg time, tn — (A/2) (c^ 2 + c^ 2 ) c t /h, where 
A is average billiard area, and h is finite-difference step 
length, was then computed by directly solving the gov- 
erning ODEs The procedure was performed for a 
number of realizations in order to obtain ensemble av- 
erage of the spectral energy. Expected linear growth in 
time of the energy peaks centered at combinatoric fre- 
quencies, wi±2 = |wi±W2|, was observed. By varying 
central frequencies of the source peaks, and calculating 
average energy growth rate of the combinatoric peaks, 
the coupling strength 10 was recovered 

N(wi ±2 ,wi,W2) = [1/2) E x ± 2 ulul/D{u x ± 2 )E x E 2 , (5) 

here E\. 2 and E\±% are total energies carried by the 
source and combinatoric peaks, respectively. Theoreti- 
cal estimate of N was also obtained by means of eq. 10 
for the purpose of comparison with the DNS. Nonlinear 
operator containing specifics of nonlinearity model was 
taken in accordance with the one used in eq. (@J): 

A^a— {x,i}/3— {x' ,fc}7— {x" ,m} 

= £^i 3 kimnS (x - x') S (x - x") d 3 /dxj dx'idx" n . 
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Short-range behavior of the spatial modal correlator was 
approximated by its behavior in boundless medium, 

K°° a = { X ,i}0={ X +T,j} = A' 1 [c~ 2 + C^ 2 ] 

xUSij/2) [cf Jo (for) + cr 2 J (k t r)] 
+ (5ij/2 - hfj) [c^ 2 J 2 {k ir ) - c7 t 2 .h (ktr)] J. (6) 

Area of the billiard A enters the above expression due 
to normalization of the modes to unity. The short-range 
correlator, as it is, leads to divergence of the integrals at 
large separations in eq. (pj . To account for finite domain 
size and to keep integrals convergent spatial correlation 
of the modes was restricted by billiard diameter L: K — 
K°° exp (— r/L). The coupling strength was then found 
to scale asymptotically as \L with linear system size. 
The scaling is an artifact of problem dimensionality, and 
is not found in 3D, where N remains O jfj. With 
non-trivial dependence of the coupling strength on the 
system size, specifics of the average domain shape become 
important. They could systematically be accounted for, 
but such undertaking lies beyond the scope of the present 
work [16] . 

Comparison of the DNS results with theoretical es- 
timates is given in Fig. Energy growth rate 
at the sum frequency was utilized to obtain numer- 
ical values of the coupling strength; calculations at 
the difference frequency were carried out as well to 
check proposed symmetry of N. To factor out depen- 
dence of the coupling strength on global parameters, 
and isolate behavior associated with frequency interre- 
lation of the coupled modes, its normalized version was 
deemed best suited for analysis: N (ttii/wi+z, fci+2-k) = 

ir- 2 e- 2 A 2 [l + ( Q /c t r 2 ] 3 (fcifc 2 fc 1+2 r 4/3 N, with k = 

uj/c t . Plotted in this form the coupling strength reveals 
its increase at the resonance frequency ratios in the vicin- 
ity of u>i /o>i+2 =(l±Ci/ci)/2. The ratios correspond to 
collinear interaction of two constituent transverse waves 
producing a longitudinal one. They define a lower and 
upper bound of the source to target frequency ratio of 
the coupled modes for which the wavevectors of the con- 
stituent waves can sum, i.e. their interaction is possible 
in accordance with pseudo-momentum conservation law. 
Another constituent-wave interaction possible in physi- 
cal solids, with ci/ct > 1, involves a longitudinal and a 
transverse wave producing a longitudinal wave. Collinear 
wavevector summation for this interaction type occurs 
at wi/wi+2 = [1 ± (3 - ci/ct) I (1 + ci/ct)] /2. However, 
this resonance does not manifest itself for the given non- 
linearity model. Under normal conditions it is expected 
to be noticeably smaller than the one involving two trans- 
verse source waves due to equipartition of the energy car- 
ried by two different types of constituent waves inside a 
solid. Since most of the energy, in particular, a fraction 
(ci/ct) > 1, is carried by transverse waves, their partic- 
ipation ratio in the coupling strength is expected to be 
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Figure 1: Normalized coupling strength for (a) ci/ct = 2, 
loi+2 = 1.2ct/h; and (b) ci/c t = 3, uji+2 = l.Oct/h. Solid 
and dash lines represent theoretical estimates, symbols A and 
□ provide DNS data for nominal system sizes of 128 and 256 
grid points, averaged over 25 and 10 realizations, respectively. 
Error bars correspond to one standard deviation. 



higher than that of longitudinal ones. 

Although good qualitative, if not quantitative, agree- 
ment between theoretical estimates and DNS was found 
for parameters used, several remarks on their discrep- 
ancies ought to be made. A difference between back- 
ground off-resonance coupling strength levels, especially 
at low frequency ratios, is noted. The difference can be 
attributed to the fact that though theoretical predictions 
based on eq. © are fit-parameter free, some ambigu- 
ity in defining effective volume (2D area) of the billiard 
interior, where nonlinear mode coupling takes place, ex- 
ists. In the present theoretical estimates it was taken 
exactly equal to the average domain area. No account 
of boundary and confinement effects was made by the 
assumed correlator form © with qualitative long-range 
correlation dependence on the order of domain size. Near 
the boundaries increased correlation of the modes due to 
boundary conditions is found, greater coherence in their 
interaction and corresponding increase in the observed 
coupling strength is expected. Such effects can play a no- 
ticeable role for system sizes within current DNS reach. 
In addition, at low frequency ratios wavelength of the 
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Figure 2: Normalized coupling strength for ci/c t = 2, 
oji+2 = 1.2ct/h. Dot, solid and dash lines represent theo- 
retical estimates at uii/uii+2 = 0.15, 0.25 and 0.35. Symbols 
A, x and □ provide DNS data for nominal system sizes of 
128, 192, 256 and 384 grid points, averaged over 25, 25, 10 
and 10 realizations, respectively. 

lowest source mode becomes large with respect to bound- 
ary roughness, and comparable to domain size. Regu- 
lar structure of the mode acquired in this case leads to 
deviation from the assumed modal statistics, and is ex- 
pected to provide increased coherence of mode coupling 
in the interior of the billiard. Large-scale symmetry of 
the square billiard, however, seems to be of less signifi- 
cance in this DNS calculations of the coupling 
strength for ensemble of rough-boundary 3:2 aspect-ratio 
quarter-stadium billiards produced compatible results. 

As mentioned above, coupling strength is predicted to 
asymptotically grow as \[L with system size. The pre- 
diction, however, is not supported by the DNS data, see 
Fig. The origins and implications of the disagreement 



are not fully understood at the time. However, it is spec- 
ulated that it can be the result of frequency smoothing 
imposed on the coupling function by the finite width of 
the source energy distribution. The smoothing leads to 
the fact that only N, averaged over characteristic peak 
width, could be recovered with the help of eq. (0. Both 
decrease and broadening of the mode coupling resonance 
peak will follow from this smoothing. System sizes cur- 
rently available for DNS realizations do not allow direct 
observation of the expected vL scaling, or draw a con- 
clusive statement regarding its absence. Decrease of the 
source energy distribution width would be problematic, 
in that it needs to contain many modes in order for sta- 
tistical approach to be valid. 

In summary, existence of the wavevector resonance 
in the mode coupling strength, responsible for non- 
linear redistribution of the average spectral energy, 
was verified by direct numerical simulation in a two- 
wavespeed chaotic billiard with nonlinearity representa- 
tive of isotropic elastic solids. Qualitative agreement 
with theoretical estimates in the absence of fit param- 
eters was observed. Location of the coupling strength 
peak was found to correspond to collinear interaction of 
constituent waves of different wavespeeds. The location 
was determined from arguments based on conservation 
principles, and to leading order depends solely on linear 
system quantities. 
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